/**************************************************************************
** Qt Creator license header template
**   Special keywords: Programmer 2012/6/6 2012
**   Environment variables: 
**   To protect a percent sign, use '%'.
**************************************************************************/

#ifndef VOXELIZATION_H
#define VOXELIZATION_H
#include <vector>
#include <assert.h>
#include <vcg/complex/exception.h>
#include <vcg/container/simple_temporary_data.h>
#include <vcg/space/box3.h>
//#include <vcg/complex/allocate.h>
#include <vcg\complex\complex.h>
#include <vcg/complex/algorithms/update/bounding.h>
#include <vcg/complex/append.h>
#include <vcg/space/intersection3.h>

template <class MeshType>
class Voxelization
{
public:

    typedef typename MeshType::VertexType            VertexType;
    typedef typename MeshType::VertexPointer         VertexPointer;
    typedef typename MeshType::VertexIterator        VertexIterator;
    typedef typename MeshType::ConstVertexIterator   ConstVertexIterator;
    typedef typename MeshType::EdgeIterator          EdgeIterator;
    typedef typename MeshType::EdgePointer           EdgePointer;
    typedef	typename MeshType::ScalarType		     ScalarType;
    typedef typename MeshType::FaceType              FaceType;
    typedef typename MeshType::FacePointer           FacePointer;
    typedef typename MeshType::FaceIterator          FaceIterator;
    typedef typename MeshType::ConstFaceIterator     ConstFaceIterator;
    typedef typename MeshType::FaceContainer         FaceContainer;
    typedef typename MeshType::CoordType             CoordType;

private:
    MeshType mesh;

public:
    Voxelization(){}

    struct Resolution{

        unsigned int xSteps;
        unsigned int ySteps;
        unsigned int zSteps;

        ScalarType xStepLength;
        ScalarType yStepLength;
        ScalarType zStepLength;
    };

    vcg::Box3f& getBBox()
	{
        return mesh.bbox;
    }

    int getVertexNum()
	{
        return mesh.VN();
    }
	int getFaceNum()
	{
		return mesh.FN();
	}

    void initialize(MeshType& m)
	{

        vcg::tri::Append<MeshType,MeshType>::MeshCopy(mesh,m);
        VertexIterator vIterator = mesh.vert.begin();
        for(;vIterator!=mesh.vert.end();vIterator++)
            vIterator->P() = m.Tr * vIterator->P();

        vcg::tri::UpdateBounding<MeshType>::Box(mesh);
    }
    void voxelize(unsigned int res,std::vector<bool>& boolList,std::vector<vcg::Box3<ScalarType> >& voxels,bool isVoxels = false,bool isFast = true)
	{


        Resolution resBox;        
        ScalarType x = mesh.bbox.DimX();
        ScalarType y = mesh.bbox.DimY();
        ScalarType z = mesh.bbox.DimZ();

        //deal with the plane
//        x<=1e-10?(x = 0.0;resBox.xSteps = 1;):(resBox.xSteps = res;);
//        y<=1e-10?(y = 0.0;resBox.ySteps = 1;):(resBox.ySteps = res;);
//        z<=1e-10?(z = 0.0;resBox.zSteps = 1;):(resBox.zSteps = res;);
        if(x<=1e-10){x = 0.0;resBox.xSteps = 1;}else{resBox.xSteps = res;}
        if(y<=1e-10){y = 0.0;resBox.ySteps = 1;}else{resBox.ySteps = res;}
        if(z<=1e-10){z = 0.0;resBox.zSteps = 1;}else{resBox.zSteps = res;}

        resBox.xStepLength = x/(ScalarType)resBox.xSteps;
        resBox.yStepLength = y/(ScalarType)resBox.ySteps;
        resBox.zStepLength = z/(ScalarType)resBox.zSteps;


        vcg::Point3<ScalarType> trans = mesh.bbox.min;

        vcg::Box3<ScalarType> faceBox;
        for(FaceIterator fi = mesh.face.begin(); fi != mesh.face.end(); ++fi)
		{
            faceBox.SetNull();
            faceBox.Add(fi->P(0));
            faceBox.Add(fi->P(1));
            faceBox.Add(fi->P(2));

            vcg::Point3<ScalarType> min = faceBox.min;
            vcg::Point3<ScalarType> max = faceBox.max;
            unsigned int voxelIndex = 0;
            unsigned int startX = 0;
            unsigned int endX  = 0;
            unsigned int startY = 0;
            unsigned int endY = 0;
            unsigned int startZ = 0;
            unsigned int endZ = 0;

            if(resBox.xSteps == 1)
			{
                startX = 1;
                endX = 1;
            }
			else
			{
                startX = (unsigned int)((min.X()-trans.X())/resBox.xStepLength);
                endX = (unsigned int)((max.X()-trans.X())/resBox.xStepLength);
            }

            if(resBox.ySteps == 1)
			{
                startY = 1;
                endY = 1;
            }
			else
			{
                startY = (unsigned int)((min.Y()-trans.Y())/resBox.yStepLength);
                endY = (unsigned int)((max.Y()-trans.Y())/resBox.yStepLength);
            }

            if(resBox.zSteps ==1)
			{
                startZ = 1;
                endZ = 1;
            }
			else
			{
                startZ = (unsigned int)((min.Z()-trans.Z())/resBox.zStepLength);
                endZ = (unsigned int)((max.Z()-trans.Z())/resBox.zStepLength);
            }

            boolList.resize(resBox.xSteps*resBox.ySteps*resBox.zSteps,false);
            vcg::Box3<ScalarType> voxelBox;
            for(unsigned int ix = startX;ix<=endX;ix++)
			{
                for(unsigned int iy = startY;iy <= endY;iy++)
				{
                    for(unsigned int iz = startZ; iz <= endZ; iz++)
					{
                        voxelIndex = ix*resBox.zSteps*resBox.ySteps + iy*resBox.zSteps + iz;
                        if(voxelIndex >= (resBox.zSteps*resBox.xSteps*resBox.ySteps))
                            continue;

                        voxelBox.min = vcg::Point3<ScalarType>((float)ix*resBox.xStepLength+trans.X(),
                                                               (float)iy*resBox.yStepLength+trans.Y(),
                                                               (float)iz*resBox.zStepLength+trans.Z());

                        voxelBox.max = vcg::Point3<ScalarType>((float)(ix+1)*resBox.xStepLength+trans.X(),
                                                               (float)(iy+1)*resBox.yStepLength+trans.Y(),
                                                               (float)(iz+1)*resBox.zStepLength+trans.Z());

                        if(isFast)
						{
                            if(!boolList[voxelIndex] && faceBox.Collide(voxelBox))
                                boolList[voxelIndex] = true;
                            
                        }
                        else
						{

                            if(!boolList[voxelIndex]&&vcg::IntersectionTriangleBox<CMeshO::ScalarType>(voxelBox,fi->P(0),fi->P(1),fi->P(2))/*!boolList[voxelIndex] && faceBox.Collide(voxelBox)*/)
                                boolList[voxelIndex] = true;
                  
                        }
                        if(isVoxels)//save the voxel box for rendering.
                            voxels.push_back(voxelBox);

                    }//end iz
                }//end iy
            }//end ix

        }//end for each face


    }//end voxelize


    //if the similarity is bigger than ratio, the two model are same.
    static bool compare(std::vector<bool>& blst1,std::vector<bool>& blst2,ScalarType ratio)
	{
        bool isSame = false;
        unsigned int count = 0;
        if(blst1.size() != blst2.size())
		{
            return false;//because num of voxels is different.
        }
        for(size_t i = 0; i < blst1.size(); i++){
            if((blst1[i]&&blst2[i]) || (!blst1[i]&&!blst2[i]))
                count++;
        }

        ScalarType rat = (ScalarType)count / (ScalarType)blst1.size();

        if(rat >= ratio)
            isSame = true;

        return isSame;
    }

    static bool isSameBBox(vcg::Box3<ScalarType>& bbox1,vcg::Box3<ScalarType>& bbox2)
	{
        ScalarType dist1 = vcg::Distance<ScalarType>(bbox1.max,bbox2.max);
        ScalarType dist2 = vcg::Distance<ScalarType>(bbox1.min,bbox2.min);

        if(dist1 <1e-2 && dist2 < 1e-2)
            return true;
        return false;
    }

};
#endif // VOXELIZATION_H
